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SUMMARY 

A computational  model  has  been  developed  for  the  turl  ulerit 
wake  left  ly  a 1 ody  moving  through  a stat  ly  stra.tifieri  medium. 
Details  of  the  wake  growth,  collapse  and  generation  of  internal 
waves  were  examined  ly  the  application  of  a second-order  closure, 
approach  to  turlulent  flow  developed  at  A.R.A.P.  o :er  the  past 
few  years.  Predictions  of  the  model  have  heen  verified  ly 
comparison  with  a wide  variety  of  wake  flows  including  wakes 
with  no  mcmentufm,  wakes  with  axial  momentum,  wakes  with  angular 
momentum,  and  for  wakes  in  both  stratified  and  unstratified  fluids. 

The  influence  of  ambient  density  gradient,  initial  density 
perturbation,  axial  momentum,  propeller-induced  swirl,  and 
vertical  lift  forces  are  all  investigated.  A numl er  of  model 
runs  demonstrate  that  the  primary  variable  affecting  the  strength 
of  the  generated  internal  waves  is  the  initial  Richardson  numl  er, 
with  the  first  local  maximum  of  the  vertical  height  of  the  wake 
scaling  inversely  with  the  l/8th  power  of  the  initial  Richardson 
number.  The  decay  rate  of  the  turbulence  intensity  as  a function 
of  distance  normalized  with  respect  to  the  l ody  diamete  appears 
less  sensitive  to  the  Froude  number  than  it  is  to  .momentum* 


Part  II  of  this  report  reviews  the  numerical  procedures 
involved  in  the  actual  computational  scheme,  and  includes  a 
complete  listing  of  the  VIA KE  program,  plus  a summary  of  its 
one rat ion. 
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NOMENCLATURE 


a,A,b,c 

C. 
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r 
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Re 
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model  constants 

drag  coefficient  of  generating  body 

porosity  coefficient  in  eq. 

diameter  of  generating  body 

QE  isotropic  correction  factor 

Nonzero  forcing  function  of  the  Poisson 
eq.  (2.  IE)  for  7. 

waxe  Frouae  number  - U/rnN 

gravitational  acceleration 
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Brunt -Vaisala  frequency  - [ — ( S , ) d,  Q oz  ] s 

pressure 

Prana tl  number  - k/v 

square  root  of  twice  the  turbulent  Kinetic 
e ne  r gy  normalized  by  u 

wake  radius 

initial  wake  radius 
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its  maximum  vaxue 

Reynolds  number1  = Ur.  v 

j i1 

9 r r' 

Rlcharason  number  of  turbulence  = r7Nc 

i max 
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affects  solution 

model  constant 
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S1,S2,S3'S4 
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ui,uJ,uk 

U 

v 

vc 

w 

WD 

Xi,XJ,Xk 
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constants  of  v.he  scale  eq.  (2.21) 
t line 

velocity  departure  in  free  stream  directlon 
normallzed  by  U 

Cartesian  velocity  components 

free  stream  uniform  velocity 

horizontal  velocity  normalized  by  U 

model  constant 

vertical  velocity  normalized  by  U 

streaming  velocity  defect  at  the  centerline 

Cartesian  coordinates 

coordinate  in  free  stream  direction 

coordinate  in  horizontal,  vertical  direction 
normalized  by  initial  wake  radius 
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Dirac  delta  function 

coefficient  of  laminar  diffusion  for  density 
perturbat ion 

microscale  of  turbulent  dissipation  = 

A/ ( a + bq.'./v  ) 2 

macroscale  of  turbulent  model 
macroscale  in  the  y,  z directions 
kinematic  viscosity 

_ /-\ 

perturbation  pressure  = [ p + j gpQdz  ] ,-y^U^ 
density,  normalized  density 

normalized  perturbation  density  = ( P ~P0 ) /( r^ dp  /dz 
ambient  i'luid  density 
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superscripts 
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stream  function  - / vdz 
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denotes  time  average 

denotes  fluctuation  about  the  mean  value 


denotes  maximum  value 
denotes  initial  conditions 
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1.  INTRODUCTION 

The  passage  of  a tody  through  a density  stratified  medium 
creates  a fluid  dynamical  problem  of  great  complexity.  The 
tody  motion  generates  a turbulent  wake  which  is  capable  of 
transforming  turbulent  kinetic  energy  to  potential  energy  by 
interacting  with  its  surrounding  density  environment.  The 
potential  energy  is  in  turn  eventually  transformed  into  internal 
gravity  waves  which  may  radiate  away  from  the  immediate  surround- 
ings. The  numerical  simulation  of  this  wake  collapse  Is  of 
special  concern  in  this  final  report. 

Other  investigators  have  studied  the  collapsing  wake  from 
experiments  (ref.  1-3)  and  simple  theoretical  models  (refs.  4-11), 
The  approach  at  A.R.A.P.  has  been  to  work  with  the  full  equations 
of  motion  for  an  incompressible  Boussinesq  fluid  and  apply  the 
technique  of.  invariant  second-order-closure  (ref.  12)  to  the  i 
dynamic  turbulence  equations.  This  approach  simulates  the 
turbulent  wake  w.y  three  mean  velocities,  a perturbation  density, 
perturbation  pressure,  and  ten  turbulence  correlations. 

Since  the  full  scope  of  the  problem  is  admittedly  immense, 
we  have  chosen  to  move  slowly  into  the  full  solution  for  the 
fifteen  unknowns.  Our  first  work  involved  verification  of  the 
invariant  technique  for  turbulent  flow  in  a density  stratified 
medium  \y  applying  it  to  the  constant  shear  stress  sublayer  in 
the  atmosphere  (ref.  13).  This  work  was  followed  by  a numerical 
simulation  of  the  axisymmetric  flow  of  a jet,  the  momentumless 
wake  of  Naudasc'ner,  and  the  wake  of  Chevray  containing  signifi- 
cant mean  momentum  in  the  streaming  direction  (ref.  14)  . We 
then  began  the  development  of  a three-dimensional  steady  flow 
computer  program  that  would  enable  us  to  simulate  nonaxisymme trie 
flows.  First  results  from  that  program  in  application  to  flow  /’ 
regimes  restricted  to  the  presence  of  either  a perturbation  in  1 
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streaming  velocity  or  the  cross-plane  velocities  (but  not  both) 
were  reported  soon  afterward  (ref.  15).  This  report  details  the 
final  restructuring  of  the  computer  program  to  handle  all  three 
momentum  equations  simultaneously,  and  discusses  the  application 
of  the  full  equation  set  to  laboratory  data  and  Initial  condition 
variability.  In  Section  2 we  review  the  governing  equations, 
analyze  the  treatment  of  a dynamic  scale  equation,  and  discuss 
the  numerical  techniques  employed  in  the  solution  of  the  derived 
equations.  In  Section  3 ’;t  make  comparisons  of  laboratory 
observations  for  both  stratified  and  unstratified  wakes  with 
A.R.A.P.  numerical  simulations.  In  Section  4 we  study  the 
sensitivity  of  various  flow  parameters  to  changes  in  initial 
conditions . 
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formulation  of  the  model 


:?.a  Governing  Equations 

For  an  3 ncompressible  boussinesq  fluid  in  the  presence  of 
turbulence,  all  flow  parameters  may  be  written  as  the  sum  of 
mean  and  fluctuating  parts.  The  time -averaged  equations  of 
motion  become 


du*  u ! a . du.  v 

i j o ( i \ 

+ 377  l v sr:  ; 

x)  J J 


(2.2) 


When  the  background  ambient  density  gradient  is  constant* 
the-  diffusion  equation  for  the  perturbation  density  K*  = p - , ^ 
becomes 


2.3) 


Subtracting  the  mean  motions,  using  eqs.  ( 2. 1 ) -( 2. 3 ) , from 
the  full  differential  equations,  we  may  then  derive  exact  equa- 
tions for  the  Reynolds  stress  correlation  u • u | end  the  correla- 
tions  involving  the  density  fluctuation  , 1 . These  become 
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To  use  Eq . >,2.3)  for  a variable  ambient  density  gradient  it  would  be 
necessary  to  assume  that  d/dx,  k 0,^/dx^  were  negligible. 
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The  closure  of  eqs . (2.1)-(2.5)  is  obtained  by  deriving 
(or  assuming)  relationships  between  the  unknown  third-order 
correlations  and  the  second-order  correlations  and  mean  flow 
gradients.  The  procedure  involved  and  the  models  currently  used 
are  discussed  in  greater  detail  in  refs.  12  through  15.  At  this 
Juncture  it  is  sufficient  to  present  the  modeled  equations  and 
the  mean  flow  equations  in  normalized  form.  Hereafter,  lengths 
are  normalized  by  the  initial  wake  radius  t velocities  by 

the  free  stream  velocity  U in  the  marching  direction  x ; and 
perturbation  density  by  -r^d^/'dz  . In  this  way  we  arrive  at 
three  nondimens ional  parameters:  Reynolds  number.  Re  = Ur./v; 

O O _ "I  / r-  -J- 

Froude  number,  Fr  = [ ( -gr^  dp  /dz)/p  IT ) / ^ ; and  Pranatl 

number,  Pr  = </v . The  pressure  occurring  in  ec.  (2.2)  is  replaced 
by  a normalized  perturbation  pressure  ( p + P grQdz  I/PqU^.  By 

assuming  that  velocities  in  the  wake  are  close  *io  the  free  stream 
velocity;  and  neglecting  terms  of  order  (wv/U)^,  the  continuity 
equation  reduces  to  the  two-dimensional,  cross-plane  continuity 
equation.  This  may  be  used  In  combination  with  the  divergence 
of  the  cross-plane  momentum  equations  to  generate  a Poisson 
equation  for  the  perturbation  pressure.  In  eq.  (2.4)  we  may 

set  1 = J and  perform  the  summation  to  obtain  the  turbulent 

2 

energy  q = u 'u ' + v 1 v ' + w 'w 1 The  principal  equations  then 
become,  with  - d'Vdy^  + d^/dz^: 
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The  modeled  form  of  the  complementary  Reynolds  stress  equations 
become 
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The  microscale  X in  eqs.  (2.7),  (2.13),  (2.14)  and  (2.15)  is 
related  to  the  macroscale  A through  the  expression 


X = A /(a  -i-  bq  A/v ) 


The  turbulence  parameters  a,  A,  b,  s,  and  vc  as  determined 
in  refs.  13  and  14  have  values  of  2.5-  0.75?  0.125?  1.8  and  C. 
respectively.  Before  discussing  the  approximations  needed  to 
solve  cqs.  (2.7}-(  2.15)  numerically  on  A.R.A  P.'s  computer 
system,  we  first  turn  to  a discussion  of  the  ietermi.iat ion  of 
the  scale  length  A. 


2,b  Treatment  of  the  Turbulent  Macroscale 

In  all  the  modeling  worK  that  A.R.A.P.  nad  done  up  to  1974, 
we  assumed  that  the  behavior  of  the  macroscale  /.  depended 
only  upon  the  gross  features  of  the  particular  problem  being 
addressed.  As  the  flow  regime  sniftea  (say,  from  axisymmetric 
flow  to  flat  plate  boundary  layer  flow),  the  rules  governing 
'•'he  determination  of  /.  changed,  but  a consistent  pattern  of 
requiring  that  A be  proportional  to  the  spread  of  the  mean 
velocity  or  turbulent  distribution  qc  was  maintained.  For  the 
axisymmetric  flow  configurations,  a good  estimate  for  A was 
obtained  by  setting 

* 

A - cr 


(2.16) 
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£ 

where  r defines  the  distance  from  the  center  r = 0 of  the 

p 

flow  co  the  radius  where  q has  dropped  to  a quarter  of  its 
maximum  value.  The  constant  c has  a value  of  0.2  (ref.  14). 

Thus  A is  set  equal  to  the  constant  value  given  by  eq.  (2.16) 
across  the  entire  flow  width. 

Our  analysis  of  the  Monin-Obukuov  sublayer  (ref.  13)  indicates 
that  in  a stratified  fluid,  the  local  scale  - while  being  governed 
overall  by  eq.  (2.16)  - is  also  restricted  in  its  maximum  value 
by  the  critical  Richardson  number  Ri*  of  the  flow,  so  that 


A < 

/ Ri  q 

X1Z 
\ lo  6Z 

It  < 0 

(2.17) 

o 

where  q~  and 

are  found 

locally . 

Equation  (2.17) 

then  limits 

the  value  of  A 

in  a region 

of  stable 

stratification; 

the  current 

value  of  Ri*  = 0.25. 

For  two-dimensional  cross-plane  motion  in  the  y and  z directions 

within  a stratified  turbulent  wake,  it  seems  logical  to  examine  the 
o 

q decay  in  both  directions  and  apply  eq.  (2.16)  to  determine 
and  Az  . Since  the  wake  collapses  in  an  elliptic -like  behavior, 
we  have  chosen  to  write  the  appearing  in  the  equation  for  \ 

and  the  Isotropy  terms  in  the  turbulence  equations  as 


(2.1d) 


2 

where  A is  first  obtained  from  the  q decay,  then  restricted 
if  necessary  by  eq.  (2.17). 


Even  though  the  gross 
with  experiments  (witness 
desire  to  investigate  the 


scale  A has  given  us  reasonable  agreement 
refs.  12-15) , there  still  remains  the 
dynamic  behavior  of  the  scale  based  on 
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the  solution  of  its  own  differential  equation*. 


Expressions  for  a dynamic  scale  equation  have  been  formulated 
by  several  observers  starting  from  the  two-point  velocity  correla- 
tions (refs.  12,  16  and  17).  Others  have  begun  with  the  vorticity 
fluctuations  or  the  dissipation  function  (refs.  l8  and  19).  As 
Bradshaw  (ref.  20)  and  Mellor  and  Herring  (ref.  21)  point  out,  the 
principal  terms  in  the  resulting  A equations  are  essentially 
the  same  for  all  of  these  formulations,  namely 


DA 

Dt 


A vw,  ^ 

S-,  U»U1  -r— — - S0v  — % + DIFFUSION  TERMS 

* q J J A 


(2.19) 


where  the  constant  s^  multiplies  the  production  term  and  s^ 
multiplies  che  dissipation.  The  principal  difference  in  expressions 
lies  in  the  construction  of  the  turbulent  diffusion  terms.  After  a 
great  deal  of  numerical  investigation,  we  feel  that  an  adequate 
dynamic  scale  equation  should  contain  at  least  two  diffusion  terms 

DIFFUSION  TERMS  - |-  (,/.  ) - »4  f 

1 i i 

where  the  first  term  comes  from  our  traditional  definition  of 
turbulent  diffusion,  found  in  eq.  (2.7)  ana  following  equations, 
and  s^  multiplies  a term  which  permits  the  diffusion  of 
to  depend  on  the  derivatives  of  q . More  complicated  expressions 
may  be  picked  for  the  diffusion,  buc  only  at  the  expense  of 
introducing  more  constants  to  be  determined.  For  stratified  flows, 
a term  proportional  to  (g^/,o)  u^t  1 should  also  be  added  to 
complete  the  current  selection  for  the  scale  equation.  Assembled, 
these  terms  give  the  following  equation  for  the  normalized  dynamic 
scale  A 


* 

The  analysis  of  a dynamic  scale  equation  was  funded  jointly  at 
A.R.A.P.  by  N00019-72-C-0413  and  the  Environmental  Protection 
Agency  under  EPA  68-02-1310. 
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To  obtain  a consistant  scale  equation  applicable  to  many  differ- 
ent flow  definitions,  several  different  well-defined  experiments 
must  be  matched  concurrently  if  the  resulting  coefficients  are 
to  have  any  invariant  validity. 


The  coefficient  s^  may  be  estimatea  from  the  dee^y  nf 
homogeneous  grid  turbulence.  If  homogeneous  turbulence  is 
assumed  to  decay  as 
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-n 
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in  the  limit  of  Re 


then  the  q equation,  eq.  (2.7),  gives 


A 


so  that  eq.  (2.21)  gives 

3 o 


U dA  _ n-2 
bq  dx  n 


(2.23) 


(2.24) 


A recent  review  of  grid  turbulence  experiments  by  Gad-el-Hak. 
and  Cori'Sln  (ref.  22)  shows  values  of  n predominately  between 
1.0  and  1.3  with  many  values  lying  near  1,?5.  This  value  of  n 
gives  a value  of  s^  = - 0.6. 

Relationships  between  the  other  three  scale  constants  may 
be  found  by  investigating  flows  near  a boi  idary  (ref.  13).  In 
steady,  neutral  flow  near  a wall  we  know  that 


A = oz 


(2.25) 
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where  z is  the  distance  normal  from  the  wall  and  a Is  a 
proportionality  constant  currently  set  at  O.65.  In  this  constant 
shear-stress  layer,  eq.  (2.21)  coupled  with  the  energy  equation 
yields  the  relationship 

s4  = vc  + \ (8i  “ s2^  (2.26) 

In  the  case  of  a stably  stratified  layer,  where  the  Richardson 
number  is  equal  to  its  critical  value  and  both  A and  q approach 
constant  values,  analysis  demonstrates  that  for  eq.  (2.21)  to  be 
consistent  with  eq.  (2.17)  requires  that 

s3  = S1  + ^S1  ” 82^  (2.27) 

Equations  (2.26)  and  (2.27)  give  us  two  relationships  between 
Sp  s2,  s0  and  s^.  Numerical  computer  fit  with  axisymmetric  Jets 
and  wakes  (with  or  without  momentum  ) , plus  a flat  plate  boundary 
layer,  and  a neutral  planetary  boundary  layer,  have  been  used  to 
provide  estimates  of  the  one  remaining  constant  and  give  a compatible 
middle  ground  across  which  the  dynamic  scale  equation  can  reasonably 
govern.  Our  best  fit  to  date  gives  s^  = - O.35,  s^  = 0.8  and 
= 0.375. 

Before  we  detail  the  comparison  of  the  dynamic  scale  with 
several  experimental  flows,  it  is  necessary  to  examine  the  appli- 
cation of  appropriate  boundary  conditions  to  eq.  (2.21).  Although 
in  most  flow  problems  the  boundary  conditions  for  the  mean  vari- 

p i 

ables  such  as  or  u^  are  known,  the  proper  or  accurate  boundary 
condition  on  A may  be  indeterminate . In  our  numerical  investi- 
gation of  eq.  (2.21),  the  edge  condition  on  \ was  taken  as 
proportional  to  the  size  of  the  gross  scale  found  by  the  straight- 
forward techniques  discussed  earlier.  Thus,  for  example!  we 
expect  the  largest  eddies  to  spread  to  the  edge  of  the  axisymmetric 
flows,  and  therefore  require  that  the  edge  value  on  \ be  twice 
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the  gross  scale  found  by  eq.  (2.16).  A similar  application  is 
made  to  the  other  flows  we  have  considered,  with  different 
factors  multiplying  the  gross  scale. 

i 

Our  comparison  of  model  predictions  for  axisyinmetrlc  flows 
made  in  ref.  14  may  now  be  recalculated  using  eq,  (2.21)  rather 
than  the  gross  constraint  eq.  (2.16).  For  the  free  jet  calcula- 
tion the  Initial  condition  on  A is  unimportant  since  only  the 
final  self-similar  results  are  compared  with  observations.  That 
comparison  is  made  in  figs.  2.1  and  2.2  with  the  results  of 
ref.  14  and  the  observations  of  V/ygnanski  and  Fiedler  (ref.  2'?9. 
Figure  2.3  compares  the  self-similar  radial  distributions  of  A 
with  the  longitudinal  integral  scale  measurements  of  Wygnanski 
and  Fiedler.  The  distributions  ?r>.»  seen  to  be  quite  similar. 

For  the  wake  observations  of  Chevray  (ref.  24)  and  Naudaseher 
(ref,  25)  we  assume  that  \ is  determined  initially  by  eq.  (2.16). 
The  decay  of  some  of  the  Important  wake  parameters  are  shown  in 
figs.  2.4  through  2.7.  Use  of  the  dynamic  scale  equation  does 
not  materially  affect  the  results.  The  radial  distributions 
obtained  some  distance  downstream  of  the  initial  station  are 
given  in  fig.  2.8.  The  dashed  line  reflects  the  approximation 
used  previously.  Relatively  little  difference  is  observed  between 
the  two  calculations. 

A comparison  of  the  effect  of  the  dynamic  scale  equation  in 
the  flat  plate  boundary  layer-  with  the  experimental  data  compiled 
by  Coles  (ref.  26)  is  shown  in  figs.  2.9  and  2.10.  Also  shown  on 
this  figure  are  the  results  from  the  gross  scale  assumption  of 
eq.  (2.16).  The  flow  profiles  possess  very  little  difference. 

We  may  also  compare  our  model  predictions  for  entrainment 
rates  with  the  experiment  of  Kato  and  Phillips  (ref.  27).  They 
measured  the  entrainment  of  water  in  an  annular  tank  with  a shear- 
ing stress  applied  at  the  water's  surface  by  a rectangular  mesh 


Figure  2.1:  Comparison  of  self-similar  jet  velocity  defect 

and  shear  stress  data  ofWygnanski  and  Fiedler 
(ref.  23)  fcith  model  predictions  for  scale  as 
determined  by  eq.  (2.16)  ( — ) or  by  eq.  (2.21) 


SCALE  VARIATION  FOR  A FREE  JET 


Figure  2.3:  Comparison  of  dynamic  scale  variation  with  radius 

for  the  self-similar  Jet.  Model  prediction  per 

eq.  (2.16)  ( — ) and  eq.  (2.21)  ( ):Wygnanskl 

and  Fiedler  (ref.  23)  longitudinal  integral  scale 
variation,  o . 


Comparison  of  the  downstream  decay  of  velocity  defect, 
shear,  and  turbulent  energy  with  scale  from  eq.  (2.16) 

( — ) and  from  eq.  (2.21)  ( ).  Initial  conditions  and 

comparison  data  points  as  given  by  Chevray  (ref.  24). 


Naudascher  , x/OsIOO 
Chevray,  x/Ds68 


Figure  2.8:  Comparison  of  scale  behavior  across  the  turbulent 

v;ake  as  given  by  eq.  (2.16)  ( ) and  by  eq.  (2.21) 

( — , ) for  the  conditions  of  Naudascher 

(ref.  2h)  and  Chevray  (ref,  24). 
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Figure  2.10:  Comparison  of  scale  behavior  in  predicting  the  law  of  the  wall 
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flyscreei.  The  rate  of  spreading  of  the  entrained  fluid  into  the 
linear  iensity  field  of  the  salt  water  could  be  measured  visually 
by  the  Injected  dye  position.  Their  experimental  results  of 
entrainment  rate  vs.  flow  Richardson  number  are  given  in  fig.  2.11. 
Our  comparisons  using  eq.  (2.16)  and  eq.  (2.21)  are  given  by  the 
solid  and  dashed  curves,  respectively.  Although  the  dynamic  scale 
rounds  the  apparently  straight-line  prediction  of  the  gross  scale 
model,  it  does  not  materially  change  the  behavior  of  the  entrainment 
with  flow  Richardson  number.  Neither  model  predicts  quite  as  steep 
a decrease  of  entrainment  with  Richardson  number  as  the  data  would 
indicate,  but  we  judge  either  prediction  to  be  satisfactory. 

In  summary  then,  it  appears  as  though  eq.  (2.21),  constrained 
by  boundary  conditions  tied  to  the  gross  spread  of  the  turbulence, 
provides  a reasonable  estimate  for  A.  However,  the  behavior  of 
A does  not  depart  strongly  from  that  by  the  gross  scale  features. 
Most  of  the  calculations  made  in  this  report  were  made  with  A 

p 

tied  to  the  gross  features  of  the  q profile  and  the  Richardson 
cutoff.  The  dynamic  scale  was  used  in  the  investigation  of  the 
sensitivity  of  the  wake  to  the  initial  scale  presented  in  Section  4. 
The  implementation  of  the  dynamic  scale  equation  in  a previously 
unexamined  flow  configuration  must  be  done  carefully. 

2.c  Computational  Approximations 

Equations  (2.7)  - (2.15)  present  a formidable  solution  task. 

To  keep  the  problem  within  a manageable  range  for  our  in-house 
computer  facility,  we  have  made  a number  of  assumptions  to  simplify 
the  numerical  analysis  while  hopefully  retaining  the  physics  of  the 
wtke  collapse. 

The  first  is  the  Quasi-Equilibrium  (QE)  approximation  (discussed 
more  fully  in  ref.  15)  uf.ed  in  our  stratified  wake  runs.  Here  we 
assume  that  the  convective  and  diffusive  character  of  all  turbulence 
properties  are  represented  by  the  dynamic  behavior  of  the  turbulent 
kinetic  energy,  eq.  (2.7).  Thus,  for  purposes  of  determining  the 
relationships  between  the  second-order  correlations,  we  discard  the 
total  derivatives  on  the  left  sides  of  eqs.  (2.1?)  - (2.15),  and 
the  diffusive  terms  on  the  right  sides.  For  the  large  Reynolds 
number  flows  that  we  consider,  \ may  be  carried  to  the  limit  of 


Figure  2.11:  Comparison  of  prediction  of  entrainment  rate  E 

vs.  shear  layer  Richardson  number  Ri0  for  scale 

as  given  by  eq.  (2.16)  ( — ) and  by  eq.  (2,21)  (— 
Data  as  given  in  ref.  27. 
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In  order  to  carry  a separate  equation  for 
algebraic  relationships  for  u ' u ' , v 1 v 


as  well  as  the 
, we  include 


and  w'w' 

in  the  energy  component 


an  isotropic  correction  factor  f(x,y,z) 
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equations  so  that  q will  always  result  from  summing 
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The  final  form  of  this  algebraic  set  of  ten  equations  becomes 


0 = - 


uiuk 


P- 


UJUh 


ujp' 
Fr  J 


Fr  * 


3.  (Tfrnj 

■ i j-  J 


3 


')  “ 2(b-f)  | ^ 


'U 


(2.29) 


(2.30) 


o 


(2.31) 


It  car.  be  seen  from  eqs.  (2.29)  - (2.31)  that  the  turbulent 

fluctuations  depend  only  upon  mean  flow  derivatives  and  the  local 
2 

q value.  It  would  seem  reasonable  to  upgrade  the  solution  scheme 

— 5 

by  retaining  the  full  dynamic  equation  for  t ' as  suggested  by 
Mellor  and  Yamada's  (ref.  28)  expansion  about  isotropic  turbulence. 
This  has  not  been  done,  however.  Solution  comparisons  with 
eas.  (2.29)  - (2.31)  in  operation,  as  opposed  to  the  use  of  the 
-ull  equations,  were  presented  for  the  axisymmetric  Naud&scher 
case  in  ref.  15.  The  agreement  in  the  dynamic  response  is  good, 
with  a maximum  difference  of  10%  occurring  in  the  mean  velocity 
departure  when  it  has  decayed  an  order  of  magnitude  from  the 
initial  conditions.  Tne  dramatic  reduction  in  the  number  of 
differential  equations  requiring  full  numerical  solution  appears 
to  give  ample  justification  for  this  slight  loss  in  accuracy. 

The  algebraic  solution  of  eqs.  (2,29)  - (2.31)  for  the 
correlations  as  a function  of  q , A , Fr  , and  the  mean  flow 
derivatives  is  not  a trivial  problem,  however.  To  further 
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facilitate  our  obtaining  a solution,  we  retain  only  the  principal 

A 

derivatives  of  u and  p in  y and  z . This  assumption 
requires  that  v and  w derivatives  be  « q/A  . However,  in 
order  not  to  lose  the  major  Influence  of  v and  w derivatives 
on  the  correlation,  the  contributions  of  'dv/dz  and  bw/dy 

for  isotropic  turbulence  are  added  to  that  obtained  by  the  above 
procedure.  The  expressions  for  the  correlations  as  used  in  our 
QE  approximation  are  given  in  Part  II. 

As  noted  in  ref.  15  we  have  also  experimented  with  modifying 
the  diffusion  terms  in  eq.  (2.7)  to  allow  turbulent  diffusion  to 
be  more  anisotropic.  This  is  done  by  replacing  v qA  by 

V 

3v  , (v  ' v '/q)  A for  diffusion  in  the  y direction  and  3v  ( w 'w  '/q 
c y c z 

in  the  z direction,  with  A , A computed  as  given  in  Section  2.b. 

y 2 

The  slight  difference  in  entrainment  this  yields  for  the  stratified 
experiment  of  Kato  and  Phillips  (ref.  27)  is  shown  in  fig.  2.12. 

The  stratified  watte  runs  reported  in  Section  A were  made  using  the 
anisotropic  form. 

Equations  (2.7)  - (2.1^)  and  the  scale  eq.  (2.21)  are  solved 
by  the  ADI  method  of  Placeman , Rachford  and  Douglas(refs  29  and 
30).  A more  complete  discussion  of  our  application  of  this 
technique  is  given  in  Part  II  of  this  final  report.  In  summary, 
the  equations  are  written  In  finite  difference  form  and  solved 
cn  an  unequally  spaced  y - z mesh,  with  the  solution  marching 

p 

downstream  in  x . A Phase  I calculation  involves  only  a , 

A 

p and  u (no  collapsing  cross -velocities ) . Phase  II,  reported 
on  In  ref.  15,  solves  for  q^  , , v , w and  v ; while  Phase  III 

encompasses  the  complete  solution  to  all  equations.  The  equations 
are  linearized  about  their  previous  x step  in  such  a way  that 
each  equation  is  solved  Implicitly  without  coupling  to  the  other 
equat ions . 

The  solution  is  relatively  slow  (on  the  order  of  6 minutes 
of  Digital  Scientific  Corporation  META-4  computer  time  per  step) 
because  of  the  need  in  the  cross-plane  for  solution  of  the 
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Figure  2.12:  Comparison  of  entrainment  rate  prediction  for  isotropic 

( — ) and  anisotropic  ( ) diffusion.  Data  (o)  from 

Kato  and  Phillips  (ref.  27). 
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pressure  perturbation  v , To  solve  the  Poisson  equation  V v - F 

at  every  step  in  the  flow,  we  add  two  additional  terms  to 

eq.  (2.12).  We  first  realize  that  our  numerical  scheme  cannot 

in  itself  maintain  a divergence-free  cross-plane  flow;  numerical 

roundoif  error  prevents  the  exact  solution  of  eq.  (2.1)  for  v 

and  w . To  aid  in  maintaining  accuracy,  since  a divergent  velocity 

will  stimulate  an  instability,  we  employ  the  standard  procedure 

(ref.  29)  of  adding  a term  of  the  form  d/dx  (7  • v)  to  the  right 

side  of  eq.  (2.12).  At  every  step  we  correct  for  zero  divergence 

at  the  next  step.  In  order  to  solve  Poisson's  equation  in  the  elliptic- 

like  region  of  a submarine  wake,  we  add  a time-like  term  to  the  left 

side  of  eq.  (2.12)  to  arrive  at  a diffusion-like  parabolic  equation 

for  r.  of  the  form 

If  = ^ - F *!x  (If  ♦&)  <2-32> 

The  pressure  perturbation  is  obtained  by  iteration  toward 
the  next  step  a Ax  distance  away,  using  the  present  pressure 
profile  as  a first  guess.  We  have  modified  the  "ideal"  time 
stepping  for  a rectangular  flow  region  (ref.  30)  to  take  steps  as 

At(i)  = yszs/PNU)  (2.33) 

where  yD  and  z_  are  the  spreads  of  the  mesh  solution  in  the 

0 5 

y and  z directions,  and  P^(i)  is  an  array  of  constants 

PN(i)  - (10.0,  30.0,  60. C,  100.0,  150.0)  (2.3A) 

Certain  minimum  error  and  cutoff  criteria  may  terminate  the 
pressure  solution  prior  to  its  completing  the  optimized  five 
iterations  (ref.  31).  Step  sizes  are  limited  by  the  possibility 
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that  the  flow  changes  too  rapidly  for  the  pressure  to  maintain  a 
divergence-free  environment.  Loss  of  numerical  accuracy  and 
solution  believablllty  follow  rapidly.  Many  checks,  both  numerical 
and  visual,  accompany  the  run  of  a collapsing  wake  beyond  one 
Brunt-Vaisala  period. 

The  principal  product  of  a collapsing  wake  is  radiated  internal 

gravity  waves.  Near  the  edges  of  our  mesh,  we  may  generally  impose 

2 A 

zero  boundary  conditions  on  q ana  u . To  require  that  , , t?  , 

v and  w also  equal  zero  here  requires  that  we  construct  an 
artificial  liner  to  contain  and  absorb  any  waves  reaching  it.  For 
this  purpose  we  add  a decay  term  -kv  to  the  right  size  of  eq.  (2.10) 
and  -kw  to  the  right  side  of  eq.  (2.11),  where 

f 0 r < sp 

k = j 5 [exp  ( +T ) + exp  (-T)]  -1  (2.33) 

- T = cp(r2  - sp2)  r > sp 

A figure  illustrating  the  successful  use  of  this  liner  is  contained 
in  ref.  15.  The  liner  cannot  lie  so  close  as  to  dictate  the  solu- 
tion, nor  can  it  be  so  far  away  that  the  mesh  loses  a proper  defini- 
tion of  the  turbulent  wake  dynamics.  As  long  as  the  liner  is 
placed  more  than  10A  away  from  the  center  of  the  wake,  the 
position  of  the  liner  alters  the  flow  solution  by  no  more  than  5$ 
of  local  maximum  values  (ref.  13). 

Taken  together,  the  numerical  simplifications  do  not  materially 
detract  from  the  ability  -„f  the  wake  program  to  perform  to  expecta- 
tion. In  the  next  Section  we  will  catalogue  some  of  the  further 
comparisons  with  experimental  data  of  our  ax i symmetric  programs  and 
the  wake  program  Itself.  The  Section  following  that  will  discuss 
the  sensitivity  studies  we  have  carried  out  with  our  operational 
computer  programs. 
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3.  COMPARISON  OF  MODEL  PREDICTIONS  WITH  LABORATORY  DATA 
3, a Unstratified  Flows 

Having  addressed  ourselves  to  the  derivation  of  the  equations 
of  motion  governing  a developing  wake  in  a stratified  fluid,  we 
now  turn  to  further  verification  of  our  model  predictions  with 
data  obtained  in  the  laboratory.  We  have  previously  (in  refs.  13 
and  14)  dealt  with  our  model  predictions  in  the  atmospheric  surface 
layer  and  the  axlsymmetric  flow  configurations  of  a free  jet  and 
wakes  with  varying  momentum.  We  will  first  investigate  some 
further  comparisons  with  unstratified  flow  experiments,  and  then 
turn  to  the  stratified  towing  tank  experiments  of  Flow  Research, 

Inc . 

Recent  wind  tunnel  work  by  Schetz  and  Jakubowski  at  VPI  (ref.  32) 
reexamined  the  N'audascher  experiments  with  a nonswirling  injection 
model  and  a swirl-induced  propeller  model  The  data  given  to  us 
was  recorded  at  their  first  station  behind  the  body,  at  x./D  = 2. 

In  the  nonswirling  case,  the  data  exhibited  some  irregularity  and 
required  interpolation  near  r = 0 to  yield  an  adequate  set  of 
Initial  conditions  with  which  to  begin  an  axlsymmetric  numerical 
calculation  as  shown  in  fig.  3*1.  Comparisons  between  the  data 
and  our  wake  predictions  are  presented  in  fig3.  3.2  - 3.4  for 
the  centerline  velocity  excess  and  the  maximum  values  of  the  shear- 
ing stress  and  longitudinal  turbulent  correlation.  No  changes 
have  been  made  In  the  model  constants  or  solution  technique  from 
that  presented  previously. 

The  centerline  ''eloclty  gives  a good  representation  to  the 

data.  It  should  be  noted  that  the  initial  data  station  at  x/D  = 2 

is  not  plotted  in  fig.  3.2  because  w^  < 0 there.  The  shear 

stress  In  fig.  3.3  matches  the  data  points  quite  well.  Figure  3,4 

plots  the  comparison  of  our  prediction  of  w 1 w 1 with  the  data 

max 

taken  In  the  experiment.  Here  we  are  at  variance  with  the  experi- 
ment. Our  theoretical  prediction  comes  relatively  quickly  to  a 
decay  rale  similar  to  Naudascher's  experimental  result  We  have 
no  explanation  for  why  the  turbulence  should  decay  so  much  slower 
in  this  momentumless  experiment. 
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Figure  3.1:  The  initial  computer  profiles  for  velocity  defect,  shear,  and  turbulent 

energy  for  the  nonswirling  data  of  Schetz  and  Jakubowskl  (ref.  32). 


Figure  3.2:  Comparison  of  the  prediction  of  the  decay  of  the 

velocity  defect  for  the  initial  conditions  given 
in  fig.  3.1  with  the  observations  (o)  of  Schetz 
and  Jakubowskl  (ref.  32). 
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Figure  3.3:  Comparison  of  the  prediction  of  the  decay  of  shear 

(stress  for  the  initial  conditions  given  in  fig.  3.1 

with  the  observations  (o)  of  Schetz  and  Jakutowskl 
(ref.  32). 
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Figure  3 . ^ : Comparison  of  the  prediction  of  the  decay  of  turbulent 

energy  for  the  initial  conditions  given  in  fig.  3.1 
with  the  observations  of  Schetz  and  Jakubowsk.1  (ref.  32). 
Also  shown  is  the  decay  characteristic  from  the  data  of 
Naudascher  (ref.  25). 
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Schetz  and  Jakubowsk.3  also  studied  a swirling  case,  and  our 
comparisons  with  their  results  are  shown  in  fig.  3.5  - 3.3. 

To  obtain  our  predictlo  '3,  we  have  used  the  A.R.A.P.  3'B  vortex  pro- 
gram for  axisymmetr 1c  flow  developed  by  R.  D.  Sullivan  for  ARL 
(ref.  33).  Our  fit  to  the  conditions  at  x/D  = 2 are  given  in 
fig.  3.5.  A dynamic  scale  equation  is  not  used  here,  but  the 
gross  scale  features  are  prescribed  as  in  the  axisymmetric,  non- 
swirling  flows . Beginning  with  the  limited  data  from  Schetz 
(rather  than  both  u' v ' and  uTwT  only  the  principal  shear  stress 
was  reported)  we  predict  the  solutions  for  centerline  w^  , 
u 'w ' and  w ’w ' „ as  shown  by  solid  lines  ir  figs.  3-6  - 3.8. 

In  fig.  3.6  we  see  that  w^  follows  the  data  fairly  well.  The 
comparisons  with  uTwT  and  w 1 w 1 are  not  as  favorable.  The 

shear  stress  is  predicted  to  decay  faster  than  that  observed 
while  the  axial  velocity  fluctuation  decays  slower  than  observed. 

At  least  the  asymptotic  rates  of  decay  of  the  axial  velocity 
fluctuations  are  the  same. 

Overall,  the  comparisons  with  Schetz ’s  data  are  somewhat 
disappointing.  We  expect  the  model  to  be  more  accurate  than 
these  comparisons  Indicate.  Whether  these  discrepancies  Indicate 
pome  defect  in  the  model,  such  as  the  treatment  of  the  scale  A 
(particularly  In  the  swirling  case)  or  some  inconsistency  in  our 
interpretation  of  the  data,  cannot  be  determlnea  until  we  have 
the  opportunity  to  examine  a more  complete  report lng  of  the 
experimental  observations. 

While  at  TRW,  Gran  (ref.  3;4)  performed  a reexamination  of 
Naudascher's  experiments  on  momentumlesR  wakes,  but  with  the 
presence  of  swirl.  Our  comparisons  with  chese  results  are 
presented  in  figs.  3.9  - 3.12.  In  fig.  3.9  we  do  quite  well 
in  our  estimate  of  centerline  velocity  excess;  we  are  just  below 
the  error  bars  near  x/L  = 10  . Figure  3. 10  compares  our  predic- 
tion of  the  maximum  axial  turbulent  intensity;  we  asymptote 
properly  but  lie  slightly  above  Gran’s  data.  Our  prediction  of 
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Ax’ol  position,  x/D 

Figure  3.9:  Comparison  of  the  prediction  of  enc  ’ecay  of  centerline  velocity  excess 

with  the  data  of  Gran  (ref.  3*0  and  Naudascher  (ref.  25). 


Figure  3.11:  Comparison  of  the  prediction  of  the  decay  of  the  swirling  velocity  with 

the  data  of  Gran  (ref.  34). 


Comparison  of  the  prediction  of  the  turbulent  wake  radius  decay  with 
the  data  of  Gran  (ref.  3*0  and  Naudascher  (ref.  25). 
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the  maximum  swirling  velocity  in  fig.  3.11  is  quite  good;  our 
prediction  for  the  wake  radius  - in  fig.  3.12  - lands  well  Inside 
the  data  points.  Taken  together,  this  comparison  with  swirl  is 
much  better  than  our  comparison  with  the  swirl  data  of  Schetz. 

Recent  data  from  Flow  Research,  Inc.  (ref.  35)  dealing  with 
the  Fr  =-  10.2  stratified  wake  experimental  run,  also  included 
an  examination  of  the  unstx\atified  experiments  of  Chevray.  Our 
model  predictions  for  this  case  are  given  in  figs,  3.13  - 3.16. 

Our  comparisons  with  maximum  velocity  defect  w^  , maximum 

vertical  and  horizontal  turbulent  intensity,  and  the  mean 
radius  r*  are  all  good.  We  conclude  from  these  very 
favorable  comparisons  in  unstratified  flows  that  our  turbulence 
model  and  gross  scale  equation  are  giving  us  a good  picture  of 
the  actual  physics  taking  place  within  these  flow  configurations. 
This  assurance  permits  us  to  investigate  problems  for  which  simple 
comparisons  are  not  possible;  certainly  much  of  our  stratified 
flow  research  falls  into  this  category. 

3.b  Stratified  Flows 

Quite  favorable  comparisons  between  model  predictions  and 
Wu's  strong  collapse  data  for  a fully  mixed  wake  (ref.  3)  and 
with  Hartman  and  Lewis'  linear  analysis  of  a collapsing  wake 
(ref.  5)  have  been  made  previously  (ref.  15).  We  include  here 
comparisons  with  some  of  Flow  Research,  Inc.'s  data.  This  has 
been  done  in  two  ways;  by  making  one  particular  run  corresponding 
to  initial  wake  conditions  supplied  to  us,  and  by  comparing  our 
vertical  scale  behavior  with  the  behavior  of  several  FRI  flow 
visualization  experiments  (ref.  36).  Since  adequate  data  to 
compare  model  predictions  with  experimental  data  for  the  one 
particular  run  has  not  been  supplied  to  us,  the  initial  condi- 
tions and  the  model  predictions  are  presented  in  the  Appendix. 

A comparison  between  model  predictions  of  vertical  scale 
behavior  and  the  behavior  of  many  FRI  flow  visualization  experi- 
ments is  shown  in  fig.  3.17.  The  solid  lines  are  an  accumulation 


Figure  3.13:  Comparison  of  the  prediction  of  the  decay  of  velocity  defect  with  the 

drag-body  experiments  of  FRI  (ref.  35). 


Figure  3.15:  Comparison  of  the  prediction  of  the  streaming  turbulent  energy 

component  with  the  FRI  data  (ref.  35). 
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of  several  A.R.A.P.  computer  runs,  with  the  data  coming  from 

five  FRI  experimental  runs.  The  scale  for  the  data  points  Is  on 

the  right  while  that  for  the  model  predictions  Is  on  the  left. 

Since  we  do  not  Know  the  exact  proportionality  between  the  height 

as  defined  by  the  edge  of  the  turbulence  and  the  height  as  defined 

by  the  dye  visualisation,  nor  the  proportionality  constant  between 

Cn  and  Initial  q2„  , an  arbitrary  adjustment  between  the  two 

v max 

scales  has  been  made  to  allow  the  first  local  maximum  in  H to 

coincide  for  both  curves.  The  qualitative  agreement  between  the 

predictions  and  the  observations  is  apparent.  Tiie  lower  Ri^ 

('or  higher  Fr  ) le,  the  longer  the  wake  fellows  approximately 
1/4 

the  t / growth  law  before  reaching  its  local  maximum  at 
approximately  the  same  normalized  time.  This  curve  will  be 
studied  more  In  the  next  Section  of  this  report. 

Overall,  our  agreement  with  stratified  and  unstratifiec 
laboratory  fiov/s  gives  confidence  to  our  solution  predictions 
in  flows  for  which  verification  will  be  difficult.  We  now  turn 
to  many  of  these  types  of  flows,  as  we  study  the  sensitivity  of 
the  wake  dynamics  to  changes  in  the  wake  initial  conditions. 
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4.  SENSITIVITY  OF  WAKE  DEVELOPMENT  TO  INITIAL  CONDITIONS 

One  of  the  principal  purposes  of  the  present  report  is  a 
recording  of  the  sensitivity  of  wake  collapse  due  to  changes  in 
either  the  initial  wake  conditions  or  the  ambient  fluid.  We 
will  investigate  changes  in  the  Initial  turbulence  level  and  macro- 
scale, variations  in  the  initial  density  profile,  changes  in  axial, 
vertical  and  angular  momentum,  and  variations  in  the  ambient  density 
gradient . 

4. a Sensitivity  to  Initial  Richardson  Number 

As  seen  In  eqs,  (2.7)  to  (2.15)  there  are  only  three  dimension- 
less parameters  in  the  governing  equations,  Fr,  Re  and  Pr.  and  the 
last  two  of  these  disappear  as  long  as  qA./v  » a/b  - 20  according 
to  our  turbulent  model.  The  primary  variable  is  then  Fr.  It  is 
enlightening  to  consider  eqs.  (2.7)  to  (2.15)  with  qQ  used  to 
normalize  velocities  rather  than  1)  . The  dimensionless  parameter 

- O 

Fr  is  then  replaced  by  the  Richardson  number  R1Q  . Further,  if 
the  streamwlse  direction  x is  normalized  by  U/N  then  q /U 
also  disappears  from  the  steady  wake  equations,  i.e.,  Ud/dx  trans- 
forms to  (Hi  )1//2  d/dx . It  is  thus  clear  that  the  two  parameters 
o 

q /U  and  Fr  may  be  conveniently  combined  into  the  single 
parameter  Ri  . An  ini •; ,i.al  indication  of  the  sensitivity  of  wake 
development  to  Ri  was  presented  in  ref.  15.  Herein  we  will  add 

O >~T 

four  wake  runs  at  Ri  values  of  Ri  = 2.18  x 10  , 2.18  x 10 

o o 

0.00925  and  0.872.  We  will  discuss  their  behavior  and  present 
contoi 'S  of  important  flow  quantities  at  half  E.V.  intervals. 

It  should  first  be  stated  that  almost  alO  runs  presented  here 

p 

began  six  body  diameters  behind  the  body  with  q - 0,0108  and 

Mid  A 

u ^0.08.  These  conditions  are  approximate  1 v those  observed  in 
the  experiment  by  Gran  (ref.  34).  Unless  otherwise  stated  the 
turbulence,  and  velocity  profiles  correspond  to  the  Naudascher  profiles. 
The  initial  density  profile  is  either  zero  everywhere  (generally  v;hen 
v and  w are  nonzero)  or  equal  to  the  contrived  density  profi3.e 


; r < 1 
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i-  L.  E.  ~ (^-i) 

U exp[-2(r2-l)l ; r > 1 

* 2 
where  r = r = 1 is  the  position  at  which  q reaches  one- 

fourth  of  its  maximum  value.  This  ’’linear-exponential"  profile 
coirespcnds  to  a well-mixed  initial  wake  (since  the  oackground 
dpo/dz  = - 1 under  cur  normalization,  for  r < 1 , the  perturba- 
tion density  balances  tne  background). 

/N 

Figure  4.1  gives  the  time  history  of  p-  „ and  q „ for 
^ iriax  max 

RiQ  = 2.18  x 10  . The  v and  w maximums  are  below  scale. 

With  Hi  this  small,  the  flow  must  run  to  rather  large  x/E 
before  collapsing.  Consequently,  v and  w never  reach  large 
values,  Qmax  decays  as  the  - 3/4  power  across  the  B.V.  scale, 
and  pmax  increases  in  response  to  the  wake  spreading,  then 
decreases  near  O.'j  B.V.  and  begins  to  oscillate.  Because  "collapse" 
takes  so  long  in  physical  distance  downstream,  we  were  able  to  use 

O 

a Phase  I (q~  , p , u)  calculation  to  0.01  E.V.  before  including 
the  two  cross-plane  velocities  v and  w , or  the  necessary 
iteration  for  the  perturbation  pressure  -•  . Figures  4.2  to  4.4 
present  the  contour  pictures  of  q^  , p arid  y for  this  waae 
at  0.5,  1.0,  and  1.5  B.V.  after  generation.  The  qc  contours 
exhibit  the  spreading  of  the  turbulent  wake  in  the  y direction, 
and  its  corresponding  shrinking  in  the  z direction,  to  produce 
an  ell ipti.c -1  ike  turbulent  wake  region.  The  density  profiles 
reflect  the  mixing  of  the  heavier  and  lighter  fluids  leading  to 
collapse,  over-collapse,  and  oscillation  about  a neutral  condi- 
tion. The  streamline  patterns  for  v in  fig.  4.4  shows  the 
typical  streamline  generation  for  a collapsing  wake:  the  definite 

presence  of  one  vortex  at.  0.5  B.V.;  two  vortices  at  1 B.V.;  and 
three  vortices  at  1.3  B.V.  In  the  fringe  areas  the  collapsing  wake 
is  generating  the  next  vortex  to  enter  the  body  of  flow.  Note 
that  most  of  the  internal  wave  dynamics  is  outside  the  main  body 
of  turbulence,  yet  the-  scales  and  associated  with  she 

turbulence  are  adequate  to  define  the  gross  scales  for  the  entire 
collapse  region. 


i 


j 

! 


Figui^e  4.1: 


The  decay  of  maximum  values  of  q and  p with 
normalized  time  (one  B.V.  ^occurs  at  Nt  = 2tt) 


for 

u , 
max 


R1o  “ 


2.18  x 10-7.  q2 


* qmax  0 ' °1°8 ' t'max  - lf 
= 0.080  (Phase  I only)  and  v = w = 0.  Time 


is  measured  from  the  point  of  wake  initialization. 


I 


F 

| 
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If  we  now  increase  RiQ  to  Ri0  -=  2.18  x lO-'5,  we  obtain  the 
wake  characteristics  shown  in  figs.  4.5  - 4.8.  Note  for  fig.  4.5 
that  q again  decays  with  little  regard  for  internal  waves; 

^ uld  A 

c'max  does  seme  readjusting  before  building  and  then  falling 

near  0.5  B.V.;  and  v and  w are  present  in  our  scale 

max  max 

only  after  ~ 0.2  B.V.  The  vertical  velocity  drops  more  rapidly 
than  v as  additional  wave  modes  are  generated  by  the  collaps- 
ing wake.  The  contours  for  q^  , ~ and  ■;  all  exhibit  the  same 
general  shapes  as  for  the  smaller  Richardson  number  flow  in 
figs.  4.2-  4.4. 

If  we  now  consider  an  intermediate  value  of  RiQ  = 0.00925, 
we  obtain  the  characteristics  shown  in  figs.  4,9  - 4.13.  For 
comparison  purposes  later  on,  we  have  started  the  Phase  III  calcu- 
lation with  u velocity  and  with  7 = 0 everywhere.  The  time 
history  plot  in  fig.  4.9  again  reflects  the  general  behavior  found 
in  a wake  collapse.  The  cross-velocities  v and  v:  are  generally 

noise  until  near  0.05  E.V. 

If  we  now  jump  to  the  other  end  of  the  scale  and  compute 
the  dynamics  for  a fast  (RiQ  = 0.872)  collapse,  we  find  the 
results  plotted  in  figs.  4.14  and  4,15.  From  the  maximum  values 
In  fig.  4.14,  we  see  that  qrr.0  slowly  turns  to  begin  its  power- 
law  oecay  (present  through  the  collapses  at  the  three  smaller 
Rl_  values)  while  , is  effectively  stable  and  tails  off  to 

begin  oscillating  near  0.5  B.V.  Because  collapse  is  so  imminent, 
we  usea  a Phase  II  calculation  throughout.  The  rr.aximuns  of  v 
and  w appear  to  grow  as  the  first  power  until  the  collapse 
begins.  A cross-over  in  maximum  value  occurs,  with  w under- 

fTla  X 

going  collapse  effects  in  line  with  . while  v barely 

ind  x ma  x 

feels  the  effect.  The  contour  plots  esent  much  the  same 
picture  as  for-  the  previous  three  .uQ  runs.  We  only  show  here 
the  plots  at  1 B.V.  in  fig.  4.15. 


Maximum  valuts 
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Figure  4.14: 


The  decay  of  maximvn 
with  normalized  time 


values  of 
for  Hl0  = 


q , p , v and  w 
0.872  and 
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Taken  together,  we  see  a fairly  consistent  picture  of  wake 

2 ~ 

development:  q decaying,  p holding  fairly  constant  before 

falling;  and  v and  w growing  with  their  maximums  occurring 
near  0.5  B.V.  The  effect  of  Ri0  can  perhaps  be  best  indicated 
by  assembling  the  pl^ts  of  the  vertical  wake  height  as  a function 
of  time  for  a number  of  different  values  of  R1  . Since  the  free 

0 i 

stream  direction  x may  appropriately  be  scaled  with  RiQ2  as 

discussed  earlier  and  since  the  momentumless  wake  radius  grows 

i 

approximately  as  xu  (ref.  14)  it  is  appropriate  to  multiply 
the  vertical  height  H by  Rio1/7^  to  correlate  the  different 

runs.  Such  a plot  is  shown  in  fig.  4,l6  with  H measured  oy 

2 

the  point  at  which  q falls  to  1/4  its  maximum  value.  For 

the  smallest  Ri  (=  F.l8  x 10"^),  we  see  that  Az  grows  with 

the  l/4lh  power  of  time  or  distance,  reaching  a broad  maximum 

near  0.3  B.V.  For  larger  values  of  Ri  , the  curves  enter  almost 

o 

horizontally,  and  then  turn  to  attempt  the  quarter-power  law 
before  reaching  a maximum  at  approximately  the  same  normalised 
time.  For-  the  largest  values  of  R1q  , the  growth  phase  of  wake 
development  is  completely  missing  with  the  reduction  in  height 
coming  earlier  and  being  more  pronounced. 


In  fig.  4.16  time  is  measured  from  the  time  of  wake  initial- 
ization. In  order  to  relate  this  to  time  measured  from  generation 
it  would  be  necessary  to  add  an  Incremental  time  N£L  = (x/r^) 

(q  /U)(Ri  )5  which  for  assumed  conditions  of  r,  ^ 1-/1 


‘max/ 


ana 


‘max 


7U  - 0.104  at  6 di  ameters  behind  the  body  would  give 


Nitt  = 1,25(R1  )2  . The  horizontal  portion  of  the  curves  at 
early  times  is  caused  by  the  coordinate  stretching  Inherent  in 
the  normalization  ana  the  logarithmic  scale. 


The  qualitative  behavior  ex:  ibitea  in  fig.  4.1o  agrees  very 
well  with  experimental  observations  as  seen  in  fig.  3-17. 
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4.b  Sensitivity  to  Initial  Density  Distribution 

Figures  4.14  and  4,15  were  begun  with  our  linear-exponential 

form  of  the  initial  density  profile,  eq.  (4.1).  We  may  also 

investigate  the  sensitivity  of  wake  development  to  cnanging  this 

initial  density  condition.  If  we  begin  with  no  perturbation 

density  anywhere  with  R1  = 0.872,  we  obtain  the  decay  character- 

° o '•> 

istics  shown  in  fig.  4,17  and  the  contour  plots  for  q , p and 

^ given  in  fig.  4,18.  It  can  be  seen  by  comparison  with  the 

maximum  values  of  fig.  4.14  that  the  qmov  profile  has  hardly 

changed.  This  would  suggest  that  the  initial  density  does  not 

change  the  q^  decay,  even  through  collapse.  Of  course,  , 

begins  from  zero  value,  and  grows  at  a rate  comparable  to  the 

A 

first  power.  But  p never  gets  to  as  large  a value  as  beginning 
with  nenzero  7 due  to  the  relatively  short  time  available  prior 

A 

to  collapse.  Since  the  magnitude  of  drives  the  cross-plane 

velocities,  vmax  and  wmax  050  not  reach  as  lar£e  values  as 
before,  AIbo  note  that  the  maximum  value  of  w_,  occurs  at  a 

IllOA 

slightly  later  time. 

The  effect  of  initial  density  profile  on  the  subsequent  wake 
development  is  reduced  as  RiQ  decreases.  For  values  of  R1Q  < 0.01 
there  is  no  discernible  effect  of  initial  density  profile  on  wake 
geometry.  In  the  other  limit  as  Rio  is  increased  (corresponding 
to  decreasing  Fr  for-  a specified  body),  the  initial  conditions 
increasingly  dominate  the  wake.  Below  some  threshold  value  of  Fr 
we  expect  any  internal  waves  generated  by  the  wake  t o be  dominated 
by  waves  generated  directly  by  the  body  so  the  behavior  of  the 
wake  in  this  regime  is  relatively  unimportant  for  fur-field  calcu- 
lations. To  demonstrate  what  may  happen  on  the  border  of  such 
a low  Fr1  regime  we  have  made  a calculation  in  which  the  initial 
wake  profiles  are  dominated  by  the  body  wave  to  such  an  extent  as 
to  create  a negative  initial  density  perturbation  (i.ef,  the  wake 
has  negative  potential  energy)  equal  to  half  the  value  of  , ^ E 
across  the  profile.  These  results  are  shown  in  figs.  4.19  and  4.20. 


Maxirrrjm  values 


Nt/2  w 


Figure  4.19:  Tue  decay  of  maximum  values  of  q , 

with  normalized^ time  for  Hlc  = 0.87 

q2„  = 0.0108,  o = - 0.9,  v - w = 0 , 
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Again  we  see  that  Qmax  is  unchanged;  Pmax  is  at  a different! 

level,  but  its  decay  behavior  is  much  the  same  as  before;  while 

• ^ 

v and  w grow  and  collapse  in  a range  of  values  between  p.  = 0 

^ /N  ~ 

and  = p L.  E.  For  this  last  case  vmax  and  wtnax  in  fact 
have  slightly  different  growth  rates,  in  very  similar  manner  to 
their  growth  in  the  run  in  the  Appendix  generated  from  the  'initial 
conditions  supplied  by  the  FRI  data.  Here  vmax  grows  at  a 
slightly  higher  rate  than  wm_v  ; a cross-over  in  curves  occurs; 
and  w shows  a strong  drop  off  after  the  maximum  values  are 
reached. 

The  different  behavior  of  the  vertical  height  of  the  wake 

for  these  three  runs  are  shown  in  fig,  4.21  by  plotting  A_  versus 

normalized  time.  We  have  also  indicated  by  the  dashed  line  the 

run  made  to  compare  with  Appendix  A.  The  slight  dipping  initially 

may  mean  that  we  have  not  given  a consistent  description  of  the 

initial  profiles  to  the  wake  program  for  this  case.  The  classic 

picture  of  a collapsing  wake  (ref.  1)  is  not  evidenced- ih  any  of  these 

runs.  A maximum  vertical  height  is  reached  in  each  case,  but 

later  in  all  the  runs  we  see  A rebounding  and  increasing. 

z 

4.c  Sensitivity  to  Initial  Turbulent  Scale 

In  order  to  investigate  the  effect  of  initial  turbulent  scale 
on  wake  development,  it  is  necessary  to  use  the  dynamic  scale 
equation  discussed  in  Section  2.  Returning  to  our  axisymmetric 
program  for  the  full  set  of  equations,  we  can  test  the  consequences 
of  a variable  A across  the  wake  profile.  This  is  in  fact  the 
reason  why  we  chose  to  derive  and  study  a dynamic  scale  equation 
whose  present  form  is  given  in  eq.  (2.21). 

For  this  test  we  assume  that  A is  small  near  the  center  or 
the  watte  (to  simulate  turbulence  chopped  up  by  the  propeller)  and 
growing  linearly  out  .to  the  edge,  fig.  4.22.  A comparison  with  our 
previous  Initial  condition,  eq.  (2.16)  is  axso  given.  The  very 
s ,iall  scale  near  r = 0 causes  a sharp  reduction  in  the  turbulent 

I 

* 1 

| 


Figure  4.21:  The  behavior  of  the  vertical  scale  Az  with  normalized  time  for 

different  initial  conditions  on  p . 


jj 
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Figure  4 22:  The  initial  profile  distribution  of  scale  length 

A versus  radius  r for  an  axisymmetric  momentumless 

wake.  The  comparison  is  between  the  profile  of 

eq.  (2.16)  ( ) and  our  contrived  condition  ( ). 


correlation  there,  due  to  the  increased  dissipation  rate.  Then 

the  flow  settles  swiftly  into  its  typical  type  of  behavior. 

The  initial  small  scale  has  been  chewed  up,  and  replaced  by 

larger  eddies  associated  with  the  larger  scale  of  the  turbulence 

near  the  wake  edge.  The  maximum  values  of  w^  , w ' w ' , andu'w' 

as  a function  of  distance  behind  the  body  for  variable  and  constant 

A are  as  shown  in  fig.  4.23.  There  is  very  little  evidence  to 

suggest  that  such  a strong  difference  existed  at  iriitisttion  of'  the 

two  runs.  From  the  fact  that  the  axial  velocity  fluctuations 

agree  for  x/D  > 50  we  can  expect  that  the  initial  turbulent 

scale  will  have  no  effect  on  ^ „ for  Fr  > 40. 

max  . i 

s 

4.d  Sensitivity  to  Angular  Momentum 

We  now  turn  to  an  examination  of  the  sensitivity  of  wake 
collapse  to  variations  in  momentum  in  the  three  directions 
(marching  x ; horizontal  y ; vertical  z).  We  first  consider 
the  influence  of  an  initial  ^swirl  configuration.  The  initial 
contours  with  swirl  for  q2  , u and  ^ are  given  ir,  fig.  4.24, 

The  density  begins  with  zero  value  everywhere,  and  V'  is  the 
product  of  the  two  compatible  v and  w velocities.  The  two 
major  effects  of  the  swirl  may  be  expected  to  be  an  enhancement 
of  mixing  in  the  density  profiles  and  an  increased  radial  spread- 
ing of  the  wake.  The  first  effect  will  be  most  important  at  low 
Froude  numbers  when  the  initial  density  distribution  is  important, 
while  the  second  may  persist  even  as  Fr  -►  on  since  it  implies 
a possibly  different  decay  rate  for  the  turbulent  energy. 

Figure  4.25  shows  the  Influence  of  swirl  on  the  decay  of  the 
wake  characteristics  for  Fr  -»•  °°  when  the  initial  maximum  swirl 
velocity  is  equal  to  the  initial  velocity  defect  w^  =_0.08  U. 
Although  q decays  somewhat  slower  initially  in  -the  swirling 
case,  asymptotically  the  two  cases  nearly  parallel  each  other. 

The  ratio  of  the  swirl  velocity  to  the  velocity  defect  chosen 
for  this  sample  calculation  is  approximately  equal  to  that 
observed  by  bran  (ref.  34)  and  a factor  of  three  higher  than 
that  for  Schetz  (ref.  32)  initial  conditions.  It  thus  represents 
a relative  upper  limit  of  the  influence  that  may  be  expected  of 
pr-opelle  r- induced  swirl. 


Maximum  values 


max 


x/D 


Figure  4.25: 


The  decay  of  the  velocity  defect,  turbulent  energy 
and  maximum  swirl  velocity  for  Fr  -+  °o  for  the 

initial  conditions  of  fig.  4.24  ( ) with' that 

of  a nonswirling  wake  ( — ). 
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The  decay  characteristics  for  Fr  = 100  and  Rio  - 0.00929 

are  shown  In  figs,  4.26-4.30.  We  see  that  q^.  appears  to 

max  _o/4 

decay  slowly  at  first,  and  then  more  nearly  to  the  typical  t 
near  colxapse,  The  cross-velocities  v and  w decay  very 
similarly.  We  may  compare  f'lgs.  4.26  - 4.33  with  the  nonswirl 
run  under  Identical  conditions,  figs.  4.9  - 4.13.  We  see  that 

A 

reaches  a larger  value  ( = 1.  rather  than  0.3)  with  swirl, 

and  does  not  show  the  oscillations  present  in  the  nonswirling  case. 

The  v ana  w are  much  lower  owing  to  their  production  by 
ma x ma x 

collapse  alone,  and  not  collapse  plus  swirl. 

Decreasing  Fr  to  Fr  = 10  and  hlo  - 0.929  yields  the  plots 

in  figs,  4.31  - 4.39.  We  see  that  q , u „ , v „ and  w „ 

0 J J v nmax'  max  ’ max  max 

are  stalled  even  longer  before  they  begin  decaying.  In  fact,  it 

appears  as  though  the  collapse  dynamics  has  caught  w before  it 

decays  very  much.  Likewise,  Fmax  hoes  not  reach  as  large  a 

peak  value  prior  to  collapse.  The  contours  are  very  consistent 

here . 

We  may  also  compare  this  swirl  collapse,  fig.  4,31,  with  the 

linear-exponential  collapse  at  Ri^  _ 0 .672,  shown  earlier  in 

fig.  4.14.  Although  the  maximum  values  taste  different  routes 

to  collapse,  at  collapse  . , q „ , v „ and  w . possess 

rudA  iTicix  [iicix  rridx 

correspondingly  similar  magnitudes.  We  have  also  repeated 

fig.  4.16  here  as  fig.  4.36  with  the  two  dashed  curves  indicating 

the  two  swirling  runs. 

From  these  runs  we  conclude  that  at  low  Fr  the  major-  effect 
of  propeller-induced  swirl  appears  to  be  to  Introduce  more  complete 
mixing  Into  the  wake  as  the  wake  spreads  more  rapidly,  forcing  a 
stronger  collapse.  At  Fr  = 100  the  major  influence  appears  to 
be  an  asymmetry  remaining  in  the  flow  contours  which  permits  the 
slight  difference  in  vertical  height  seen  in  fig.  4.36  and  the 
significant  difference  in  energy  radiated  in  Table  4.1  shown  later. 
For  still  higher  values  of  Fr  we  would  expect  smaller  all  Terences 
due  to  propeller-induced  swirl. 


Maximum  values 


Nt/27T 


Figure  4.2c-:  The  decay  of  maximum  values  of  q , p , u v 

and  w for  Ri0  = 0.00925  with  swirl  and  initial 
conditions  given  in  fig.  4.24. 


'igure  4 . 2 f : 


Contours  of  constant  q' 
fig.  4.20  (see  fig.  4.2 


fo 

for 
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q2m0,  =0.0000867 
0.5  B.V. 


Figure  U.2b:  Contours  of  constant.  ( for  the  conditions  of  fig 

(see  fig.  >\.e>  for  contour  code). 


Figure  4.29:  Contours  of  constant  u for  the  conditions  oj 

Fig,  4.2o  (see  fig.  4.2  for-  contour  code). 


Maximum  values 
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Figure  4.34:  Contours  of  constant  u for  the  conditions  of 

fig.  4.31  (see  fig.  4.?  for  contour  code). 
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4.e  Sensitivity  to  Axial  Momentum 

We  may  estimate  the  effect  of  nonzero  streaming*  momentum 
(in  the  x direction)  by  returning  to  the  FRI  data  run  given 
in  the  Appendix.  If  we  use  a symmetrical  version  of  their  (\ 
and  u profiles,  along  with  our  linear-exponential  density 
profile,  we  ob tain  the  collapse  characteristics  shown  in  fir.  4.37- 
4.41.  In  comparing  the  time  history  of  the  maximum  values  with 
the  R1q  * 0.872  run,  fir.  4.14,  we  see  a very  close  similarity 
across  the  entire  E.V.  scale.  There  is  a slower  decay  in  the 
turbulence  level  reflecting  the  fact  that  a wake  with  momentum 
decays  at  a slower  rate  than  a momentunless  wake  (see  firs.  2.4  and 
2.6).  However  this  has  little  effect  on  the  cross  flow.  Appar- 
ently the  low  Fr  and  the  well-mixed  nature  of  the  initial 
condition  make  the  collapse  of  the  wake  imminent  so  the  slirht 
difference  in  turbulent  decay  rate  does  not  have  time  to  affect 
the  other  flow  quantities.  We  would  expect  the  influence  of  axial 
momentum  to  increase  as  Fr  is  increased  and  the  wake  is  permitted 
to  decay  further  before  the  major-  cross-flow  is  generated. 


4,f  Sensitivity  to  Vertical  Momentum 


When  we  include  some  lift  force,  to  simulate  vertical  momentum, 
we  begin  the  wake  run  with  profiles  in  which  we  have  two  vortices 
located  away  from  the  center  of  the  turbulent  wake  (only  the  y > 0 
side  will  be  shown).  The  maximum  velocity  of  the  lift  induced 
vortex  is  approximately  equal  to  the  axial  velocity  defect  yielding 
a positive  lift  that  is  C.5  times  the  integrated  tur't  ulent  kinetic 
energy  divided  by  the  free  stream  velocity  (for  a propulsive 
efficiency  of  0.5  this  would  correspond  to  a lift-to-drag  ratio  of 
This  ratio  is  a larger  value  than  we  would  tvnicallv  exoect 


1) 

frsrr.  a nearly  axisymmetric  iody.  The  initial  density 


V o 


zero,  and  we  run  with 


■r  =f  100,  or  Hi 


0 . 0092 


Ls  taker,  to 
The  maximum 


values  and  contour  plots  for  negative  trim  are  shown  in  fi  s.  4.42  - \ 

* 

4.-v..  We  see  that  the  maximums  exhil  it  a 1 ehavior  similar  to  the  t ehavibr 

f 


Nt/2r 


Figure  4.37: 


Decay  of  maximum  values  of  q , , u , v and 

for  R1q  = O.855  for  a wake  behind  a drag  body. 


Here,  q<L„  = 0.0108 


p = 1 

max 


u 


max 


= 0. ^05 


CO 


max 


max 


max 


- -2 

o 10  * 


H 


w 


max 


Nt/2 
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X. 


Figure 


The  decay  of  maxltrrun  values  of  q , t , u , v and  w 
for  a vortex  pair  generating  negative^,  1 if  t force  at 


= 0.00925.  Here  we  begin  with 
0 , u = 0.08  , and  lift/drag  ~ 


= 0.0108  , 


Figure  4.43;  Contours  of  constant  q2  for  the  conditions  ol 
fig.  4.42  (see  fig.  4.2  for  contour  code). 


0.5  B.V. 


r,' i cure  4 44:  Contours  of  constant  for  the  conditions  oi 

fir  4.42  (.see  fig.  4.?  for  contour  cocie)  . 


4-56 


in  the  swirling  case,  but  not  with  similar  decay  characteristics, 
particularly  for  the  three  velocity  components.  The  density  again 
shows  the  collapse  effe<_%  with  Nmax  dropping  rapidly  there. 

The  contours  exhibit  some  of  the  vortex  motion;  we  3ee  that  the 
vortex  is  eventually  dissipated  by  the  turbulence,  and  lost  by 
the  stratified  collapse. 

We  see  that  lift  force  makes  relatively  little  difference  in 
the  decay  characteristics.  If  we  now  repeat  the  calculation  but 
reverse  the  forcing  role  (to  make  the  lift  positive)  we  obtain 
the  decay  behavior  plotted  in  fig.  4.47  and  the  contours  in  figs  4.48 
4.51.  A comparison  of  fig.  4.42  and  4.47  shows  very  little  differ- 
ence across  the  B.V.  scale  for  q „ , p , u „ , and  wm_„  . 

^max  hmax  ’ max  ’ max  max 

The  contours  show  a mirror  effect,  particularly  in  u , until 
collapse  dominates  at  1 B.V. 

For  a further  clarification  of  the  role  of  lift  forces  we 
made  a run  for  the  Fr  -*•  °o  limit.  In  this  limit,  there  can  be 
no  gravity  induced  cross  flow,  but  only  the  presence  of  two  vortices 
decaying  together.  Our  calculation  shows  that  indeed  the  vortices 
drift  downward  from  their  initial  positions  and  travel  away  from 
each  other.  Such  behavior  is  well-predicted  by  linear  theory. 

In  this  case,  then,  the  density  perturbation  is  decoupled  from 
the  flow  dynamics.  In  fig.  4.52  we  show  the  printer  plot  of  density 
at  a station  x/L'  = 10.7  downstream  of  initialization.  This  type 
of  plot  Is  the  typical  output  from  the  WAKE  program  generating 
most  of  the  results  presented  in  this  report.  V.:e  can  easily  discern 
the  spirals  present  in  the  vortex  field.  Note  that  the  negative 
values  of  p have  larger  absolute  values  than  do  the  positive 
values,  reflecting  the  fact  that  the  total  wake  has  moved  downward. 

Table  4.1  presents  a summary  of  the  runs  made  investigating 
the  sensitivity  of  the  strength  of  the  gravity-induced  cross  flow 
on  initial  conditions.  In  addition  to  listing  the  figures  where 
detailed  plots  of  the  time  history  of  major  parameters  and  contours 
e.  selected  times  may  be  found,  we  have  included  two  parameters 
v'..l,n  measure  the  overall  strength  of  the  internal  waves  genex-ated 
by  ..e  wake.  These  are  the  maximum  value  of  the  cross-flow  stream 


Nt/2w 


Figure  4.47:  Decay  of  maximum  values  of  q , p , u , v and 

w for  a vortex  pair  generating  positive  lift 
force  at  Ri  = 0,00525.  Here  we  begin  with 

= 0.0108  , p = 0 , u = 0 . 08  , and  lift/drag 


u * 


Figure  4.49:  Contours  of  constant  p for  the  conditions  of 

fig.  4.47  (see  fig.  4.2  for  contour  code). 


Figure  4.50:  Contours  of  constant  u for  the  conditions  of 

fig.  4.47  (see  fig.  4.2  for  contour  code). 
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function  and  the  integrated  value  of  the  energy  radiated  away 

from  the  waRe  by  the  internal  waves.  The  first  of  these  is 

computed  directly  by  the  program  from  the  velocity  field  at  any 

time.  However  the  radiated  energy  cannot  be  computed  directly 

because  of  the  existence  of  the  absorbing  liner.  Rather,  it  Is 

computed  indirectly  from  an  energy  balance  by  subtracting  the 

energy  remaining  and  the  energy  dissipated  from  the  initial  energy. 

This  admittedly  is  not  very  accurate  when  the  energy  dissipated 

exceeds  9b/  of  the  initial  energy.  Theoretically  we  can  expect 

the  ratio  of  the  energy  radiated  to  the  initial  energy  to  decrease 

as  hi  ~ ' c as  Hi  -»  0 , while  ••  /q  r.  should  decrease  1 1 Re 
o o ’ nax  o i 


As  a further  summary,  the  curves  of  the  decay  of  the  turbulence 
intensity  as  a function  of  x,T  are  repeated  in  fig.  4.93.  The 
curves  are  labeled  as  noted  in  Table  4.1.  The  momentumless  waRes 
for  Hi  — ™ have  their  turbulence  intensity  decay  approximately 

O 

as  (x/D) "u/*4 . As  RiQ  is  increased,  the  rate  of  decay  Increases 
when  x/L'  > Fr  . When  the  curve  ends  somewhere  between  1 and  2 
Brunt -Vaisa la  periods  after  initiation,  the  turbulent  intensity 
is  approximately  20&  lower  than  the  value  of  the  KiQ  -»  <»  curve. 
Propeller  swirl  or  lift  forces  cause  some  slight  departure  from 
this  basic  trend  but  the  largest  difference  is  caused  by  the 
Introduction  of  axial  momentum.  How  the  Hi  -*  » curve  corresponds 
to  the  self-similar,  finite  momentum  decay  rate  of  (x/T)  “/  J shown 
dashed  in  fig.  4.93.  Again,  as  Fr  is  decreased,  the  decay  rate 
Of  q will  depart  from  this  when  x/L  - Fr  . This  dashed  curve 
then  represents  an  urner  bound  on  the  turbulent  intensity.  For 
waRes  which  start  almost  momentumless,  the  decay  curve  can  initially 
follow  the  steeper  momentumless  rate  and  (as  seen  in  ref.  14)  cross 
over  to  that  appropriate  for  the  momentum  waRe  far  downstream.  The 
two  momentumless  curves  shown  on  the  present  curve  for  high  Fr 
had  identically  zero  momentum  since  they  were  divided  into  a 
Phase  I and  Phase  II  calculation  as  described  earlier.  The  total 
spread  due  to  any  uncertainty  in  Fr  is  less  than  that  due  to 
uncertainties  in  momentum. 
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5.  CONCLUDING  REMARKS 

V.'e  believe  the  primary  goal  of  constructing  a practical  model 
for  wak.es  including  the  influences  of  stratification,  axial 
momentum,  ver*  leal  momentum,  and  angular  momentum  based  on  second- 
order  closure  of  the  turbulence  equations  has  been  accomplished. 
Although  the  comparisons  of  model  predictions  with  some  individual 
experimental  observations  have  been  disappointing,  agreement  with 
the  majority  of  the  experimental  data  has  been  quite  favorable. 
There  are  three  basic  areas  in  which  the  model  needs  further 
development:  1)  in  the  numerical  implementation  so  that  the  full 

dynamic  equation  for  each  Reynolds  stress  can  be  solved  when  it 
is  Indicated  that  the  quasi-equilibrium  approximation  is  not  valid; 
2)  in  further  validation  of  the  macroscale  treatment  so  that  a 
higher  confidence  can  be  placed  in  the  model  predictions  for  an 
untested  flow;  anc  3)  in  the  removal  of  the  approximation  that 
t'nc-  axial  velocity  is  a small  perturbation  from  the  free  stream 
velocity  in  the  fully  three -dimensional  wake  program  so  that 
computations  are  valid  closer  to  the  body. 

A number  of  model  runs  have  rr.-nn  made  to  investigate  the 

sensitivity  of  the  wake  development  to  Initial  conditions . 3omo 

Oi  the  more  interesting  conclusions  from  this  invostigat ion  art 

that:  1)  the  primary  variable  affecting  the  strength  of  the 

generated  internal  waves  is  the  Hie’,  ardson  number;  2)  the  decay 

of  the  turbulent  energy  as  a function  of  x/R  is  less  sensitive 

to  Froude  nu.iber  than  it  Is  to  axial  momentum;  ',)  an  Increase  in 

the  initial  density  perturbation  leads  to  an  Increase  in  the 

gravity  induced  eros3-f low  for  low  Fr  but  has  no  influence  when 

R1q  < 10  ; 4)  the  relatively  low  value  of  swirl  induced  by  a 

propeller-propelled  body  causes  a slightly  lower1  initial  decay  In 

wake  turbulence  although  the  asymptotic  decay  rate  remains  approx! - 

- ' / 4 \ 

mately  proportional  to  x J ; 5)  at  It  --  100  the  swirl  inouced 
by  a propeller  causes  an  approximate  factor  of  three  increase  in 
the  v/ave  radiated  energy  for  a wake  that  is  assumed  to  start  with 
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no  Initial  density  perturbation  primarily  oecause  of  the  more 
efficient  density  mixing  of  the  swirl  case;  6)  a lift  force 
changes  the  contour  pattern  of  the  px'imary  variables  at  1 B.V. 
without  any  significant  change  in  the  maximum  values. 
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APPENDIX 

MODEL  PREDICTIONS  FOR  STRATIFIED  V.'AKT  DEVELOPMENT 
GENERATED  FROM  FRI'S  INITIAL  CONDITIONS 


The  Initial  profiles  from  Flow  Research,  Inc.'s  stratified 
towing  tank,  experiment  (ref.  35)  were  given  at  x/D  = 6 for 

2 

several  stations  in  y and  z for  the  principal  variables  q , 

A 

p and  u . Kith  these  curves,  we  constructed  an  extrapolated 
curve  fit  which  passed  reasonably  well  through  the  giver,  profile 
a^ta.  The  initial  contours  in  the  first -quadrant  solution  plane 
are  shown  in  fig.  A.l.  An  unusual  feature  of  these  initial 
conditions  is  the  character  of  the  Initial  density.  In  our 

A 

typical  turlulent  wake  runs,  p is  a positive  departure  from 
the  background  for  z > 0 . This  is  the  behavior  expected  in  a 
wake  more  uniformly  mixed  than  the  surrounding  stratified  fluid. 
Instead,  it  seems  that  at  x/D  - 6 , the  wake  is  dominated  by 
the  body  generated  internal  wave  to  give  a wake  with  a stronger 
density  gradient  Inside  than  outside.  In  this  light,  the  FRI 
data  presents  a challenging  problem  to  the  wake  program  - to 
compute  for  a set  of  initial  conditions  its  inventors  did  not 
envision . 


Contours  of  q^  , ,,  , u , and  streamline  < are  presented 
in  figs.  A. 2 - A. 5 at  0.5,  1.0  and  1.5  brunt -Vaisala  (E.V.) 
periods  after  run  initialization.  Figure  A.o  gives  the  time 
history  cT  the  maximum  values.  Ke  observe  from  here  the  typical 
type  of  behavior:  the  decay  of  downstream;  the  spread  of 

u with  time;  and  the  buildup  of  the  cross-plane  velocities  v 
and  w to  their  maximum  value  near  0.5  b.V.  However,  because 
of  the  initial  sign  of  the-  density  perturbation,  the  development 
of  the  vertical  height  of  the  wake  is  quite  different, 
comparison  is  shown  in  fig.  c . 2 1 . 


Th  i s 


Figure  A.l: 


The  initial  contour  profiles  for  q , p and  u for  the 
Flow  Research  Inc.  data  (ref.  35)  with  Ri0  = C.855  and 

<4x  = °-009'  %ia=x-°-80  and  umax  = -305  (see  flg'  '*-2 
for  contour  code). 


Figure  A. 3:  Contours  of  constant  p for  the  conditions  of  fig.  A.l 

(see  fig.  4.2  for  contour  code). 
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Contours  of  constant  for  the  conditions  of  fig.  A.l 

(see  fig.  4.2  for  contour  code). 


Figure  A. 


imum  values 


SECURITY  CL  ASSIFICA1  ION  OF  THIS  PAOC  (Dl«l  CM*  Inll»4 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1.  REPORT  number  j.  GOVT  accession  no. 

A.R.A.P.  Report  No.  226 

3 RECIPIENT'S  CAT  ALOC  NUMBER 

« title  (tnd  tulllll,) 

Turbulent  Waives  in  a Stratified  Fluid 
Part  I:  Model  Development,  Verification, 

and  Sensitivity  to  Initial  Conditions 

s type  of  report  a period  covercd 

Final  Report 

6 PERFORMING  ORG.  REPORT  NUMBER 

> AUTHORfA) 

W.  Stephen  Lewellen 
Milton  Teske 
Coleman  duP.  Donaldson 

• CONTRACT  OR  GRANT  NUMBER(a) 

N00014-72-C-0413 

1.  PERFORMING  ORGANIZATION  NAME  ANO  AODRCSS 

Aeronautical  Research  Associates  of 
Princeton,  Inc. 

50  Washington  Rd.,  Princeton,  N.  J.  08540 

10  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  6 WORK  UNIT  NUMBERS 

DARPA  Order  No.  1910 

II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS 

Office  of  Naval  Research 

Department  of  the  Navy,  Arlington,  Va. 

12  REPORT  DATE 

August  197^ 

1 J NUMBER  OF  P AGES 
1^2 

14  MONITORING  AGInCy  name  • ADORESSCI/ dllltrtnl  from  Control. Ini  Office) 

IS  SECURITY  CLASS,  (of  1 hit  r.port) 

1S«.  DECL  ASSiEiC  ATlON'V00*NGRADlNG 
SChEOULC 

'6  DISTRIBUTION  STATEMENT  (ot  th/a  R9pott) 

■ : i -I T A 

i ...  f 

J Ap;  : i : or  public  roloaso; 

I ]'.  Ti  fni  :i  i'^d 

17.  DISTRIBUTION  ST  ATEmCnT  (oi  tha  abatract  antarad  /h  Block  20.  It  dltlarani  Itom  Report) 

'6  supplementary  NOTEi 

19  KEY  WORDS  (Com/nu*  on  tavaraa  i If  nicniiry  and  Idantlfy  by  block  number) 

Submarine  wak.es  Numerical  fluid  dynamics 

Turbulence  modeling 
Stratified  flow 
Swirling  wakes 

20  ABSTRACT  (Coniinua  on  tavaraa  alda  II  nacaaamry  and  Identity  by  block  numbar) 

A computational  model  has  been  developed  for  the  turbulent 
wake  of  a body  moving  through  a stably  stratified  fluid.  Details  of 
the  wake  growth,  collapse  and  generation  of  internal  waves  were 
examined  by  the  application  of  a second-order  closure  approach  to 
turbulent  flow  developed  at  A.R.A.P.  over  the  past  few  years. 
Predictions  of  the  model  have  been  verified  by  comparison  with  a 
wide  variety  of  wake  flows  including  wakes  with  no  momentum,  wakes 
with  axial  momentum,  wakes  with  angular  momentum,  and  for  wakes  in 

DD  I JAN^7J  1473  COITION  OF  I NOV  «S  IS  OBSOLETE 


I 


SECURITY  CLASSIFICATION  of  This  PAGE  r*>i«n  Dolt  Enlertd) 


SECURITY  CLAmnC*TlOW  0*  THI»  P«M  ■»!»»< 

both  stratified  and  unstratified  fluids.  A sensitivity 
investigation  reveals  that  the  primary  variable  affecting  the 
strength  of  the  generated  internal  waves  is  the  initial 
Pichardson  number,  with  the  first  local  maximum  of  the  vertical 
height  of  the  waVe  scaling  inversely  with  the  l/8th  power  ol 
the  initial  Richardson  number. 


jecukity  clarification  of  this  FAoenriF 


